A Jacobi algorithm for distributed model predictive control of 

dynamically coupled systems 



Dang Doan, Tamas Keviczky, Ion Necoara and Moritz Diehl 



00 

o 
o 

(N 

Oh 

in 



u 
o 



> 

cn 

d\ 
o 
oo 

o 



Abstract — In this paper we introduce an iterative Jacobi 
algoritlim for solving distributed model predictive control 
(DMPC) problems, with linear coupled dynamics and convex 
coupled constraints. The algorithm guarantees stability and 
persistent feasibility, and we provide a localized procedure for 
constructing an initial feasible solution by constraint tightening. 
Moreover, we show that the solution of the iterative process 
converges to the centralized MPC solution. The proposed 
iterative approach involves solving local optimization problems 
consisting of only few subsystems, depending on the choice 
of the designer and the sparsity of dynamical and constraint 
couplings. The gain in the overall computational load com- 
pared to the centralized problem is balanced by the increased 
communication requirements. This makes our approach more 
applicable to situations where the number of subsystems is 
large, the coupling is sparse, and local communication is 
relatively fast and cheap. A numerical example illustrates the 
effects of the local problem size, and the number of iterations 
on convergence to the centralized solution. 

I. Introduction 

Model predictive control (MPC) is the most successful 
advanced control technology implemented in industry due to 
its ability to handle complex systems with hard input and 
state constraints [5], [8], [9]. The essence of MPC is to 
determine a control profile that optimizes a cost criterion over 
a prediction window and then to apply this control profile 
until new process measurements become available. Then the 
whole procedure is repeated and feedback is incorporated by 
using the measurements to update the optimization problem 
for the next step. 

For the control problem of large-scale networked sys- 
tems, centralized MPC is considered impractical, inflexible 
and unsuitable due to information exchange requirements 
and computational aspects. The subsystems in the network 
may belong to different authorities that prevent sending all 
necessary information to one processing center Moreover, 
the optimization problem yielded by centrahzed MPC can 
be excessively large for real-time computation. In order to 
deal with these limitations, distributed MPC is proposed 
for control of such large-scale systems, by decomposing 
the overall system into small subsystems. The subsystems 
employ distinct MPC controllers, use local information from 
neighboring subsystems, and collaborate to achieve globally 
attractive solutions. 
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Approaches to distributed MPC design differ from each 
other in the problem setup. In [3], Camponogara et al. studied 
stability of coordination-based distributed MPC with several 
information exchange conditions. In [4], Dunbar and Murray 
proposed a distributed MPC scheme for problems with 
coupled cost function, utilizing predicted trajectories of the 
neighbors in each subsystem's optimization. Keviczky et al. 
proposed a distributed MPC scheme with a sufficient stability 
test for dynamically decoupled systems in [7], in which each 
subsystem optimizes also over the behaviors of its neighbors. 
Richards and How in [10] proposed a robust distributed 
MPC method for networks with coupled constraints, based 
on constraint tightening and a serial solution approach. 

A distributed MPC scheme for dynamically coupled sys- 
tems called feasible-cooperation MPC (FC-MPC) was pro- 
posed by Venkat et al. in [11], [12], based on a parallel 
synchronous approach for cooperative optimization. This 
scheme works only for input-coupled linear time-invariant 
(LTI) subsystem dynamics without state constraints, and is 
not applicable to problems with constraints between subsys- 
tems. In this paper, we propose an extension of this scheme 
in several ways in order to solve these issues. 

The distributed MPC algorithm described in this paper 
is able to handle LTI dynamics with general dynamical 
couplings, and the presence of convex coupled constraints. 
Each local controller optimizes not only for itself, but also 
for its neighbors in order to gain better overall performance. 
Global feasibility and stability are achieved, whilst the al- 
gorithm can be implemented using local communications. 
The proposed algorithm is based on an MPC framework 
with zero terminal point constraint for increased clarity and 
simplicity. We highlight an open research question that needs 
to be addressed for a full treatment of the terminal cost based 
version of this MPC framework, which would allow reduced 
conservativeness. While other distributed MPC methods typ- 
ically assume an initial feasible solution to be available, we 
incorporate a decentralized method to determine an initial 
feasible solution. 

The problem formulation is described in Section lUl fol- 
lowed by two variations of the algorithm in Section |lll] 
It is shown that an algorithm using local communication 
exists and it is equivalent to one that is based on global 
communication. In Section |IV] we analyze the feasibility, 
stability and optimality of the algorithm. Different ways of 
customizing the proposed algorithm and a trade-off between 
communications and computational aspects are discussed in 
Section |V] Finally, Section [VT] illustrates the algorithm in a 
numerical example and Section IVIII concludes the paper. 



II. Problem description 

A. Coupled subsystem model 

Consider a plant consisting of M subsystems. Each sub- 
system's dynamics is assumed to be influenced directly by 
only a small number of other subsystems. Let each subsystem 
be represented by a discrete-time, linear time-invariant model 
of the form: 



M 



^{Aijxl + Bijul), 



(1) 



where xl e 



and ul e 



are the states and control 



inputs of the i-th subsystem at time t, respectively. 

Remark 2.1 This is a very general model class for describ- 
ing dynamical coupling between subsystems and includes as 
a special case the combination of decentralized models and 
interaction models in [11]. 

We define the neighborhood of i, denoted by A/^', as 
the set of indices of subsystems that have either direct 
dynamical or convex constraint coupling with subsystem i. 
In Figure [T] we demonstrate this with an interaction map 
where each node stands for one subsystem, the dotted links 
show constraint couplings and the solid arrows represent 
dynamical couplings. The neighborhood JV^ of subsystem 
4 is the set of {4, 1, 2, 5}. We will refer to vectors and sets 
related to nodes in A/"* with a superscript +i. The collection 
of all other nodes that are not included in JV^ will be referred 
to with a superscript i. 

B. Convex coupled constraints 

Each subsystem i is assumed to have local convex coupled 
constraints involving only a small number of the others. If 
we fix the control inputs and the corresponding states of the 
nodes outside M\ the state and input constraints involving 
the nodes in can be defined in the following way: 



4' 



,M (2) 



where X~^'^(xl) and U^'''{u[) are closed and convex sets 
parameterized by the states and control inputs of nodes 
outside AA*. 

Remark 2.2 Note that the constraints involving nodes of 
A/^* in general do not depend on every other state and input 
outside A/% only on the immediate neighbors of M"^. The 
notation in Q is used for simphcity. 



C. Centralized model 



Let X — 



X 



and u = 

denote the aggregated states and inputs of the full plant, 
„i^i. vii^iw...,.^..., i=i and R-i-i=i respectively. The 
matrices A and B will denote the aggregated subsystem 
dynamics matrices and are assumed to be stabilizable: 
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Fig. 1. An interaction map showing the constraints (dotted links) and 
dynamical couplings (arrows) between subsystems. In this example, AT'* = 
{4,1,2,5}. 



The full (centralized) plant model is thus represented as: 

xt+i = Axt + But, (3) 

Remark 2.3 The centralized model defined in [3] is more 
general than the so-called composite model employed in 
[11]. In our approach, the centralized model can represent 
both couplings in states and inputs. In [11], the authors 
use an input-coupled composite model, which requires the 
subsystems' states to be decoupled, allowing only couplings 
in inputs. 

D. Centralized MFC problem 

The centralized MPC problem is formulated based on a 
typical quadratic MPC framework [8] with prediction horizon 
N, and the following quadratic cost function at time step t: 



Vt 



N-l 

E 

fc=0 



Xk,t Qxk,t + Uk,t RUk,t 



(4) 



where Xk^t denotes the centralized state vector at time t + k 
obtained by starting from the state xo,t — Xt and applying 
to system (O the input sequence uo.t, • • ■ , Ufc-i,t. Q = 
diag {Qi, • • • , Qm), R = diag (i?i, • ■ • , Km) with diag (.) 
function representing the block diagonal matrix. Matrices Qi 
are positive semidefinite and Ri are positive definite. 



Let Xt — [a 



'^w.tl 'Ut = Kt'--- .<-i,t] -The 



centralized MPC problem is then defined as: 



ultRuk,t 



(5) 



N-l 

Vt*{xt) = min V xltQxk.t 

X,,Ut -^-^ 

k=0 

s.t. Xk+i,t = Axk,t + Buk,t,k = 0, N ~1, 
Uk,t £U,k = 0,...,N ^1, 
Xk,t e X,k^l,...,N -I, 

XN.t = 0, 

xo,t = Xt, 

where U and X are defined as HiHi ^B (W"*"*) and 
nfii C-^ (-^"^O' respectively. The CE (•) operator denotes 



respectively. In other words, if C Rf then Ci? 



X Mz^i=i The vector contains the measured 

states at time step t. 

Let ~ [{uq j)-^, • ■ • , t)'^]'^ denote the optimal 

control solution of Q at time t. Then, the first sample of 
is applied to the overall system: 



(6) 



The optimization (|5]l is repeated at time t + 1, based on the 
new state xt+i- In [6] it was shown that with prediction 
horizon N long enough to allow a feasible solution to 
the optimization problem, the closed-loop system (O-© is 
stable. 

Before formulating the distributed MPC problems, we 
eliminate the state variables in the centralized MPC formu- 
lation. In the following we will also assume t = without 
loss of generality and drop subscript t for simplicity. The set 
of dynamics equations allows us to write the predicted states 
as 



(7) 



where 
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Using the above equations, we can eliminate state vari- 
ables in the centralized MPC leading to the following prob- 
lem: 

min V(u, xq) (8) 

u 

= min u^(a^Qa -t- R)u + 2(a^Q/3)^u + /3^Q/3 

U 

s.t. uGW, 

au + I3{xq) e X, 
Fm + A^xo = 0, 



where U 



n 



N 

k=l 



U and X 



n 



N 

k=l 



X. F 



[A^-^B,A^-^B, ...,B] is the last block row of a. The 
matrices Q and R are block-diagonal, built from weighting 
matrices Q and R. Therefore Q is positive semidefinite and R 
is positive definite, making the cost function ^(u, xq) strictly 
convex with any given xq. 

E. Distributed MPC problem 

We will solve problem ^ by dividing it into smaller, 
overlapping DMPC problems, with each DMPC assigned 
to one subsystem but optimizing also over neighboring 
subsystems at the same time. At each time step, by solving 
DMPCs and combining the local solutions in an iterative 
process, we will get an increasingly accurate approximate 
solution of the centralized MPC problem. 

In the DMPC problem for subsystem i, the global cost 
function is optimized with respect to a reduced set of 



variables: control inputs of i and its neighbors, denoted 
together by Each DMPC problem will guarantee that 
all constraints of the centralized MPC problem are satisfied. 
The DMPC of subsystem i can be recast into the following 
optimization problem: 

min V{u,xo) (9) 

u+» 

s.t. u eU, 

au + (3{xo) e X, 
Fu + A'^xo = 0, 



where u' " denotes the assumed inputs of all non-neighbors 
of i. For now we assume that in the beginning of each 
step, each node j transmits its assumed inputs for u-' = 
[(uq)"^, to the entire network, node i receive 

these vectors from Vj ^ A/"' to form u*'". 

Note that u is the combination of and u'. With each 
i, we can construct two pairs of matrices a' and 
F^ so that: 



au = a ' u ' - 
Fu = F+'u+' 



(10) 



By eliminating the input constraints which involves only 
u% the DMPC problem (|9]l is equivalent to the following 

min V{u,xo) (11) 

S.t. u+' e U+\ 

a+'u+' + aV + f3{xQ) G X, 
F+'u+' + F^u' + A'^xo = 0, 

in which = nf=i The optimal solution of Q, 

(fTTl l will be denoted by u+". 

For implementation, we introduce the notion of the r-step 
extended neighborhood for each subsystem i, denoted by JV^, 
which contains all nodes that can be reached from node i in 
not more than r links. A/j! is the union of subsystem indices 
in the neighborhoods of all subsystems in 7V^_i : 



(12) 



where Afl :— TV*. We see that in order to solve ( fTTT l. 
subsystem i only needs information from other subsystems 
in Nj^j^i, the initial states and assumed inputs of subsystems 
outside M^^i are redundant. 

Remark 2.4 The MPC formulation using terminal point 
constraint described above simplifies our exposition but it 
is rather conservative. This could be alleviated by using a 
dual-mode MPC formulation with a terminal cost function. 
However, in order for this to be a truly distributed approach, 
the terminal cost function associated with the terminal con- 
trollers should have a sparse structure. This would allow 
the construction of a centralized Lyapunov function in a 



local way, using only local information. In [11], the authors 
try to bypass this obstacle by using additional restrictive 
assumptions; they employ zero terminal controllers and re- 
quire all subsystems and interaction models (coupled via the 
inputs only) to be stable. These assumptions can actually 
be more conservative than using a terminal point constraint, 
preventing the application of the FC-MPC method in general 
dynamically coupled systems. Finding terminal controllers 
that lead to a structured terminal cost is an open problem 
and subject of our current research. 

III. Jacobi-type algorithm 

In this section we present an iterative procedure to approx- 
imate the centralized MPC solution by repeatedly calculating 
and combining the solutions the local DMPC problems 
described in the previous section. We will show two versions 
of our approach, which are based on Jacobi distributed opti- 
mization [2]. The proposed algorithms maintain feasibility 
of intermediate solutions and converge to the centralized 
MPC solution asymptotically. The first version uses global 
communication and can be considered as an extension of FC- 
MPC [11]. The second version relies on local communication 
and represents the main contribution of the paper We will 
show that the two versions are equivalent to each other, 
which leads to simplified analysis in Section HVl 

A. Globally and locally communicating algorithms 

For each time step t, we assume that a feasible input u-'^ is 
given for the entire system. (Section FlV-Al discusses a method 
of obtaining such a feasible initial control sequence in a 
distributed way, given a known initial condition.) In each step 
of the proposed DMPC scheme, the subsystems cooperate 
and perform a Jacobi algorithm, where each subsystem 
iteratively solves the optimization problem ( fTTT i with regards 
to its local variables, and incorporates a convex combination 
of neighboring local solutions. 

During every MPC sampling period, a distributed iterative 
loop is employed, and is indexed by p. At each iteration p, 
VL^ is updated. We will refer to vectors obtained in these 
iterations with subscript [p). 

For j3 = 1, we initialize the iteration with U(o) = u-^. Let 
u^^p^ denote the control sequence of the whole system stored 
in the memory of subsystem i at iteration p. For making 
convex combinations, each subsystem i is assigned a weight 
A* G (0, 1), satisfying X]f=i — 1- The choice of weights 
is arbitrary and could depend on the specific problem, the 
simplest choice will be equal weights (A* — We 
propose then the following iterative algorithm: 

Algorithm 3.1 (Jacobi DMPC with global communication): 
Given N, p„^ax > 0, e > 0: 

1. p ^ 1, U(o) ^ , ^ large number, V« = 1, . . . , M. 
while > e for some i and p < Pmax 
a. for each i — 1, ... ,M 

Construct new u'^'^ from U(p_i). 

Solve problem (|9]l to get uf^-* . Construct a global 



input vector u^^j from u^" and u*'". Transmit 

Ujji*^ to a central update location, 
end (for) 

b. Merge local solutions according to the following 
convex combination: 

M 

1=1 

c. Compute the progress and iterate: 

P' = -U(p-i)ll 

p ^ p + 1 
end (while) 

2. Each subsystem implements the first input value in 

3. Shift the predicted input sequence one step to make a 
feasible solution for the following MPC update: 

4. t ^ t + 1. Measure new initial states xt, go to step 1. 

Algorithm 13. 1 1 requires the existence of a central coordi- 
nator that communicates with all subsystems and performs 
the convex combination to find U(p). For implementation, 
a control scheme without global communication is desired. 
Next we introduce a variation of Algorithm 13.11 that only 
needs local communication. 

Let u* ■'^ denote the feasible input sequence of subsystem i, 
and u^l^y denote control sequence of subsystem j computed 
by subsystem i when solving its DMPC problem ( fTTT i at 
iteration p. 

Algorithm 3.2 (Jacobi DMPC with local communicatioiQ): 
Given N, pmax > 0, e > 0, and assuming each subsystem i 
knows a feasible input u-''-' for all subsystems j G Ml^j^^. 

1. p 1, «— large number, Vi = 1, . . . , M. 

while (0* > e for some i and p < Pmax 

for each i = \, . . . ,M 

a. Subsystem i solves the local problem (fTTT i. using 
{u-^'-^IVj G Nl^j^iXM"^} as assumed inputs for sub- 
systems outside TV* but inside Nj^j^i. The solution is 
comprised of {u^p", j G A/"'}. 

b. Subsystem i receives solutions for itself calculated by 
its neighbors {u'^^^*, j G A/^*}, then updates its solution 
for iterate p according to: 

"I.) - E ^Mrt + ( 1 - E "U) (15) 

c. Calculate the progress: 

d. u*'-^ <— ^\p)' subsystem i transmits new u'^-'' to all 
subsystems in 

end (for) 



p p + l 
end (while) 

2. Each subsystem i implements the first input value in 



yields u^^^ as follows: 



Uf = u, 



o,{p)- 



(16) 



3. Shift the predicted input sequence by one step to make 
a feasible solution for the following MPC update: 



^''^ = K,(p)'-" 'WAr-i,(p).0], i = l,...,M. 



4. t ^ t + 1. Measure new initial states xl, go to step 1. 

The major difference between Algorithms 13.11 and 13.21 is 
at step l.b: in Algorithm 13.11 the convex combination is 
performed on the global control input vector, while in Algo- 
rithm |32] each local controller performs convex combination 
using its local control input vectors, therefore removing the 
need of a coordinator. In the sequel, we will show that the 
two algorithms are equivalent, thus allowing us to implement 
Algorithm 13.21 while using Algorithm 13. II for analysis. 



B. Equivalence of the two algorithms 

The two crucial differences between Algorithm |3.2| and |3.1| 
are the communication requirement and the update method. 
We already mentioned that the optimization problem ( fTTT l is 
equivalent to (|9]l, thus each subsystem only has to transmit 
its new results to the subsystems inside Mjqj^i- This leads 
to the local communications in Algorithm [32] Now we will 
show that the local update of Algorithm l3.2l is also equivalent 
to the global update of Algorithm 13.11 

Consider Algorithm 13.21 at the beginning of iteration p, 
a local input vector vCp_i is given for each i. Then each 

subsystem j G TV* computes \lp^'* and sends these solutions 
to i, which forms the final update for itself Up. Note that 
i G ^ i ^ A/"*, so we have 



E - 



(P) 



E^M"(P-i)+"(P-i) (17) 



Now consider Algorithm l3.1l at the beginning of iteration 
p, starting from the same u^^ j^-j as in Algorithm 13.21 
The local problem of the j-th subsystem achieves solutions 
{u|p|i* I j G A/"*} which are equal in both algorithms and 
are transmitted to subsystem i. Then the global update U(p) 



M 



= E " 



= E k;^ + E 



M 



^ E 



E I "(P-1) 



'(p-i) 



(18) 

Comparing ( fTTI l and ( fTSl l we can see that the local update 
of Algorithm 13.21 and the global update of Algorithm 13.11 
yield the same result. This implies that Algorithm 13.21 does 
exactly what Algorithm 13 . 1 1 does, except it only needs to use 
regional information (each subsystem i needs to communi- 
cate with subsystems in the region Af^^^i). If M ^ N and 
the interaction map is relatively sparse, this region will be 
much smaller than the whole network, thus DMPC problems 
can be considered local optimization problems. 

The equivalence of the two proposed DMPC algorithms 
allows us to prove their feasibility, stability and optimality 
aspects by analyzing the globally communicating algorithm, 
which is more comprehensive than the locally communicat- 
ing algorithm. In the next sections, we refer to Algorithm l3.1l 
for analysis. 

IV. Feasibility, stability and optimality 
A. Constructing initial feasible solutions locally 

Although in current Uterature it is typically assumed that 
an initial centralized feasible solution exists and is available, 
in this section we give a simple but implementable way of 
actually constructing it in a distributed way assuming that 
the global initial state is available in advance. 

The initial feasible prediction input u-^ at time t ~ 
can be calculated locally by using an inner approximation 
of the global feasible set, which we will denote with il 
based on all the constraints appearing in dHJ and the global 
initial state, which is assumed to be available. Consider an 
inner-hyperbox (or hyperrectangular) approximation B of the 
feasible set ft, which then takes the form of a Cartesian 
product: 

B^B^ X---X B^\ 



Ben. 



(19) 



This approximation essentially decomposes and decouples 
the constraints among subsystems by performing constraint 
tightening. Each subsystem i will thus have to include B"^ 
in their local problem setup. Since the Cartesian product 
of these local constraint sets are included in the globally 
feasible set 57, any combination of local solutions within 
will be globally feasible as well. Needless to say that 



the local constraint sets that arise from this inner-hyperbox 
approximation will be in general quite conservative, but at 
the same time will allow the construction of a centralized 
feasible solution locally to initialize Algorithm 13. II 

Calculation of the inner-hyperbox approximation can be 
performed a priori and the local constraints distributed to 
each subsystem. A polynomial-time procedure to compute a 
maximum volume inner box of ft could follow the procedure 
described in [1]. Let us denote the dimension of the global 
input vector with d = J^fLi represent a box as 

B{u,u) ~ {u eM.'^ \ u < u < u}, then B{u*,u* + u*) is a 
maximum volume inner box of the full-dimensional polytope 
defined as = {w e R'' | Au < b}, where is an 

optimal solution of 

d 

max > \nvi 

n,v (20) 

s.t. Au + A'^v < b. 

The matrix is the positive part of A. Obtaining the local 
component-wise constraints is then straightforward. For 
time steps other than t = 0, we construct a feasible solution 
by performing step 3 of Algorithm 13. II 

B. Maintaining feasibility throughout the iterations 

Observe that in step l.a, we get M feasible solutions Uj-^^ 
for the centralized problem (O. In step l.b, we construct 
the new control profile U(p) as a convex combination of 
these solutions. Since problem (l8]l is a convex constrained 
QP, any convex combination of {u^^jIfZ]^ also satisfies the 
convex constraint set. Therefore U(p) is a feasible solution of 
optimization problems ([H]), and Q for all i. 

C. Stability analysis 

Showing stability of the closed-loop system (l3])-(fT4l) fol- 
lows standard arguments for the most part [6], [9]. In the 
following, we describe only the most important part for 
brevity, which considers the nonincreasing property of the 
value function. The proof in this section is closely related 
to the stability proof of the FC-MPC method in [12], the 
main difference is due to the underlying MFC schemes; this 
method uses terminal point constraint MFC while FC-MPC 
uses dual-mode MFC. 

Let pt and pt+i stand for the last iteration number of 
Algorithm 13.11 at step t and t + 1, respectively. Let Vt = 
V{u(^p^),xt) and Vt+i = y(u(p^^j), xt+i) denote the cost 
values associated with the final combined solution at step 
t and t + 1, respectively. At step t + 1, let ^Ip^i) = 

^i^lp+i)t^t) denote the global cost associated with solution 
of subsystem i at iterate p + 1, and $(p) = V{u(^p-j,xt) the 
cost corresponding to the combined solution at iterate p. 

The global cost function can be used as a Lyapunov func- 
tion, and its nonincreasing property can be shown following 
the chain: 

Vt+i<---< $(p+i) < $(p) < • • • 

• • ■ < <Vt- xfQxt - ujRut 



The two main components of the above inequality chain 
are shown in the following two subsections. 
Showing that <&(p+i) < $(p) 
The cost V{u,Xt) is a convex function of u, thus 

(M \ M 

h^E^^^«i)'^0 (21) 
1=1 / i=l 

Moreover, each is the optimizer of i-th local problem 

starting from U(p), therefore we have: 

v{»iP+iy^t)=^iP+i)^^ip)^'^'^^---'M (22) 
Substituting (l22l l into (l2Tl i leads to: 

CM \ M 

Using the above inequality, we can trace back to p ~ 1: 

Vt+i <•■•.< $(p+i) < «>(p) < • ■ • < 

Showing that <i>(i) <Vt — xjQxt ~ ujRut 
At step t + 1 and iteration p — 1, recall that the initial 
feasible solution of the centralized MFC is built by 
Algorithm 13. II at the end of step t in the following way: 

The DMFC of each subsystem i optimizes the cost with 
respect to u+* starting from u-'^, therefore Vi — 1, . . . , M: 

V(f^.xt^ <V{n^,xt) 

N-l 

fc=i 

< *(u(p,)) - xlp^Qxo,(p,) - uJ_(p^)i?uo^(p,) 
^ $(1) <Vt- xfQxt - ujRut. 

Moreover, due to the convexity of V{\i,xt) and the convex 
combination update U(i) ~ J2iLi -^''^(1*)' obtain 

/ M \ M 

*(i) = ^ E^^^'-O-^^^'^W)'"*) 

\i=l / i=l 

M M 

<^(i) < E *(i) =^ E ~ ^^Q^* - ""tRut] 
1=1 1=1 

^ *(i) <Vt~ xjQxt - ujRut 

The above inequalities show that the value function de- 
creases along closed-loop trajectories of the system. The rest 
of the proof for stability follows standard arguments found 
for instance in [6], [11]. 



D. Optimality analysis 

Using the descent approach, we will show that the solution 
of Algorithm 13.11 approaches the solution of the centralized 
MPC in (O, as p ^ oo. We characterize the optimality of the 
proposed iterative procedure by using the following results. 

Lemma 4.1: A limit point of {u(j,)} is guaranteed to exist. 

Proof: The feasible set of ^ is compact. It is shown that 
every U(p) is feasible, therefore this sequence is bounded, 
thus converges. □ 

Lemma 4.2: Every limit point of {u(p)} is an optimal 
solution of 

Proof: We will make use of the strict convexity of V{-) 
and a technique, which is inspired by the proof of Gauss- 
Seidel distributed optimization algorithms in [2]. In our 
context however, we address also the overlapping variables 
that are present in the local optimization problems. 

Let V = (v^, v^^) be a limit point of {u(p)}, assume 
{u(p')} is a subsequence of {u(p)} that converges to v. 

In the following, we drop parameter xt in V{-) for 
simplicity. Using the continuity of ) and the convergence 
of {u(p/)} to V, we see that F (u(p')) converges to V{\). 
This implies the entire sequence {F(u(p))} converges to 
V{y). It now remains to show that v minimizes V{-) over 
the feasible set of (O. 

\it]\ converges to zero. Recall 



We first show that u^l^^^ 



that at iteration p, u^^ 



(P'] 



and u 



s|l 

some e > such that llu 



forms U(p/), at iteration 
Assume the contrary, 
does not converge to zero. There exists 



and u^^^) forms uj!+,) 



s\l 

(p'+i) 



> e for all p'. 



Let us fix some 7 e (0, 1) and define 



Hp') " ^^(P') 
lies between u 



1 



and u 



(p'+l) 
(p'+ 



(23) 
only differs 



this means sj^,^ iies oeiween u(p/- 
from them in the values of u+^. 

Notice that sj^,^ belongs to a compact set and therefore 
has a limit point, denoted s^^. Since 7^1 and U(p') 7^ 
U(p!+i),Vy, we have s;^ 7^ v. 

Using convexity of V{-), we obtain 

V(sj^,)) <max{V^(u(p,)),y(uJ^,))}. (24) 



Be definition, u^^^^j^^ 
of u+^. So we have: 



V 



( s|l 



minimizes V(S) over the subspace 

<V{^\^.^<VM (25) 



From Section IIV-CI we have 

(uj+i)) < (u(p,)) , Vi = 0,---,M (26) 

M 

^(u(p.+i))<5:AT(u^i:+,)) (27) 

1=1 

Taking the Umit of ( [26] l and (|27] i. we obtain 

^lirn^ V (u(i:+i)) = F(v), V* = 0, • • • , M (28) 



As V (u(p/)) and V both converge to V^(v), 

taking limit of (|25] |. we conclude that V(y) ~ ^(sL)' for 
any 7 6 (0,1). This contradicts the strict convexity of T^(-) 
in the subspace of u+^. The contradiction establishes that 
'^(p'+i) ^'^(p') converges to zero, leading to the convergence 

of U^;,^!) to V. 

We have, by definition 

F(u(i!+i))<y(u+\u[p,)) (29) 

Taking the limit as p' tends to infinity, we obtain 

T/(v) < (u+\v^) (30) 

or V is optimizer of F(-) in the subspace of u+^. If we further 
consider V{-) in a subspace corresponding to u^, then F(v) 
is still a minimum. Thus, the necessary optimality condition 
gives VuiV^(v)'^(ui - v^) > OjVu^ £ Oi where fii is the 
feasible set of (|9]l with i = 1. 

Repeating the procedure, we obtain 



V„,F(v)^(u' -v') >0, 



(31) 



for all u* such that u* is a feasible solution of (|9]l. 

By summing up the system of equations in JsTT l for i = 
1, • • ■ , A/, we get: 

VuF(v)^(u-v) >0, (32) 

for all u that is a feasible solution of ([8]l. 

This shows that v satisfies the optimality condition of 
problem dH). □ 

Using strict convexity of it follows that v is in fact 
the global optimizer of Algorithm 13.11 

V. Communications and computational aspects 

In this section, we discuss the communications and com- 
putational aspects of our approach and illustrate the freedom 
that the designer has in choosing the appropriate trade-off 
and performance level in a certain application. 

Although the overall computational load is reduced by 
employing the distributed Algorithm 13.21 its iterative nature 
implies that communication between neighboring systems 
increases in exchange. This trade-off is illustrated in Table |l] 
which compares the communication requirements of the 
centralized and our distributed MPC approach. This overview 
suggests that our scheme is mostly applicable in situations 
where local communication is relatively fast and cheap. 

TABLE I 

Comparison of communications requirements 





Centralized MPC 


Distributed MPC 


Communication 


Global 


Local 


Each subsystem 
communicates 
with 


Central coordinator 


Other subsystems 
in {N + l)-step 
extended neighborhood 


Total number of 
messages sent 
in each time step 


2 X M 


Vmax X 2 X — 1 ^^A^-l-l 



TABLE II 

Comparison of optimization problems 



From some nonzero initial state, the system needs to be 
stabilized subject to the constraints: 





Centralized MFC 


Distributed MFC 


Number of vaiiables in 
one optimization problem 


iV X Af 


N X lA^'l 


Number of optimizations 
solved in one sampling period 


1 


Pmax X M 



< 4, = 2,...,M- 1 



(34) 



Table HI] shows the difference in size of the optimiza- 
tion problems solved by the distributed and the centralized 
method. Since jA/"*] ^ M, where M is the total number 
of subsystems, the local optimization problems in DMPC 
are much smaller than the centralized one. Note that during 
one sampling period, the local DMPC optimization problems 
are solved at most Pmax times. Nevertheless, DMPC is in 
general more computationally efficient than the centralized 
MPC, with a proper choice of Pmax- 

Choosing an appropriately high Pmax value leads to better 
performance of the whole system. The trade-off is that the 
increase of Pmax will also lead to increased communications, 
and more local optimization problems will need to be solved 
in one sampUng time. 

Another way to customize the algorithm is to expand 
the size of the neighborhood that each subsystem optimizes 
for. In the proposed Algorithms 13.11 13.21 each subsystem 
optimizes for its direct neighbors when solving the local 
optimization. We may have better performance when each 
subsystem optimizes also for its 2, 3, k-step expanded 
neighbors. Although we do not provide a formal proof of this, 
we will give an illustration in the following section through 
a numerical example. The intuition behind this phenomenon 
is nevertheless clear: each subsystem will have more precise 
predictions when it takes into account the behaviors of more 
neighboring subsystems. 

VI. Numerical example 

In this section, we illustrate the application of Algo- 
rithm 13.21 to a problem involving coupled oscillators. The 
problem setup consists of M oscillators that can move only 
along the vertical axis, and are coupled by springs that 
connect each oscillator with its two closest neighbors. An 
exogenous vertical force will be used as the control input 
for each oscillator The setup is shown in Figure |2] 

Each oscillator is considered as one subsystem. Let the 
superscript i denote the index of oscillators. The dynamics 
equation of oscillator i is then defined as 



ma' = kip'-fsv'+k2{p'^^-p')+k2{p' 



i+l 



-P 



'F\ (33) 



where p% and a' denote the position, velocity and accel- 
eration of oscillator i, respectively. The control force exerted 
at oscillator i is F' and the parameters are defined as 

ki. stiffness of vertical spring at each oscillator 

k2'- stiffness of springs that connect the oscillators 

m: mass of each oscillator 

fs'- friction coefficient of movements 



Based on dynamical couplings and constraint couplings, 
the neighborhood of each subsystem inside the chain is 
defined to contain itself and its two closest neighbors M^' = 
{i — 1, i, z + 1}, i = 2, M — 1, while for the two ends 
TV"! = {1, 2} and A/"*^ = {Af, M - 1}. We define the state 
vector as — [p'^v'^y, and the input as = F"^. The 
discretized dynamics with sampling time Tg is represented 
by the following matrices: 



— 

A-ii — 
Bn — 








Tsk2 
1 

Ts{ki-2k2 


Tsk2 




\fi = 2,...,M 

Vi = 1 



1 - TJ, 








Vj ^ i 
Vi = 1, 



The following parameters were used in the simulation 
example: 

ki = 0.4, fca = 0.3 

fs = 0.4, Ts = 0.05, m = 1 

M = 40, N = 20 

100 




= 10 



Starting from the same feasible initial state, we apply 
Algorithm 13.21 with Pmax = 2, 20 and 100. The results 
are compared to the solution obtained from the centralized 
MPC approach. The results indicate that all states of the 40 
subsystems are stabilized. Figure [3] shows the evolution of 
the overall cost achieved by DMPC compared to the cost 
of the centralized approach. We can see that the difference 
is reduced by choosing a larger Pmax value. Our analysis 
guarantees in fact that the DMPC solution converges to the 
centralized one as p tends to infinity. 

As mentioned above, another way to customize the pro- 
posed distributed MPC algorithm is for each local problem to 
consider optimizing over the inputs of subsystems in a larger 
neighborhood. Figure |4] illustrates the effect of optimizing in 
each subproblem over an r-step extended neighborhood, with 
r = {l,5,10}. Fixing the number of maximum subiterations 
to Pmax = 2, we can observe a steady improvement in perfor- 
mance until the increased neighborhood of each subsystem 
covers essentially all other subsystems and end up with a 
centralized problem. 



m 
Tf; 





















Fig. 2. Setup of coupled oscillators 



Distributed IVIPC, pmax = 2 
Distributed IVIPC, pmax = 20 

- Distributed IVIPC, pmax = 100 

- Ceritralized I\ylPC 



VII. Conclusions 

We presented a Jacobi algorithm for solving distributed 
model predictive control problems, which is able to deal 
with general linear coupled dynamics and convex coupled 
constraints. We incorporated neighboring subsystem models 
and constraints in the formulation of the local problems 
for enhanced performance. Global feasibility and stability 
were achieved, and a local implementation of the algorithm 
was given, which relies on information exchange from an 
extended set of "nearby" neighboring subsystems. It was 
shown that the distributed MPC solution converges to the 
centralized one through a localized iterative procedure. An a 
priori approximation procedure was proposed, which allows 
to construct an initial feasible solution locally by tightening 
constraints. We also discussed the trade-off between commu- 
nications and computational aspects, the effect of increasing 
the maximum number of iterations (pmax) in one sampling 
period and the potential improvements that can be gained by 
incorporating several subsystems into a local optimization. 
We are currently working on an extension of the algorithm, 
which allows the use of terminal costs in a dual-mode MPC 
formulation. 




Fig. 3. Time evolution of the global cost value of the centralized MPC in 
comparison with the distributed MPC algorithm using Pmax — 2, Prnax — 
20 and pmax = 100. 




' " Distributed MPC, r = 
- ' Distributed IVIPC, r = 

Distributed IVIPC, r = 

Centralized MPC 
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Fig. 4. Time evolution of the global cost value of distributed MPC 
algorithms with different radius of neighborhood to be optimized by one 
local controller 



